vars =  {   'Coefftrue', 'Coefficients', 'StandardErrors', 'ComputationTime',...
    'OwnElast', 'CrossElastSameSeg', 'CrossElastDiffSeg',...
    'OneOwnElast','OneCrossElastSameSeg','OneCrossElastDiffSeg',...
    'DiversionSameSeg', 'DiversionDiffSeg', 'OneDiversionSameSeg','OneDiversionDiffSeg',...
    'FunctionValue', 'ExitFlag', 'AIC', 'BIC', ...
    'SaveSegment', 'SaveX1', 'OutsideShare', 'corrDstarX1'};
load('1resultsRCLvsNLLowCorr',vars{:});
% load('2resultsRCLvsNLHighCorr',vars{:});
rc_names    = {'rc05';	'rc1';	'rc15';	'rc2';	'rc25';	'rc3';	'rc35';	'rc4';	'rc45';	'rc5'};

% I need this to save the files in excel in a loop

Excel = actxserver ('Excel.Application');
% Make sure you write the entire path
File='\\econprofiles\profiles1\n07082\Documents\MonteCarloSummaries\1Paper\1RCLvsNL\2Postanalysis\RCLvsNL2.xls';
if ~exist(File,'file')
    ExcelWorkbook = Excel.workbooks.Add;
    ExcelWorkbook.SaveAs(File,1);
    ExcelWorkbook.Close(false);
end
ExcelWorkbook = Excel.workbooks.Open(File); 

% Save results in excel

if corrDstarX1 >= 0.9
    TitleSheet = 'HighCorr';
else
    TitleSheet = 'LowCorr';
end

for rc_idx = 1:numel(rc_names);
    
    % Analysis of results
    
    % get rid of outliers (<10% only with high correlation; 0% with low correlation)
    numel(Coefficients.(rc_names{rc_idx})(:,:,:,(Coefficients.(rc_names{rc_idx})(5,1,3,:)>10)))
    Coefficients.(rc_names{rc_idx})(:,:,:,(Coefficients.(rc_names{rc_idx})(5,1,3,:)>10)) = [];

    %% Inspect results by visual look at the histogram of the coefficient distributions
    
    % Nested logit
%     figureName  = 'Distribution Pho';
%     figure2save = figure('Name',[num2str((rc_names{rc_idx})),figureName],'Color',[1 1 1],'Position',get(0,'ScreenSize'));
%     
%     histfit (squeeze((Coefficients.(rc_names{rc_idx})(4,:,2,:))));
%     saveas(figure2save,[figureName,'.bmp']);
%     normplot(squeeze((Coefficients.(rc_names{rc_idx})(4,:,2,:))));
    
    % Random Coefficient Logit
%     figureName  = 'Distribution Sigma';
%     figure2save = figure('Name',[num2str((rc_names{rc_idx})),figureName],'Color',[1 1 1],'Position',get(0,'ScreenSize'));
%     
%     histfit (squeeze((Coefficients.(rc_names{rc_idx})(5,:,3,:))));
%     saveas(figure2save,[figureName,'.bmp']);
%     normplot(squeeze((Coefficients.(rc_names{rc_idx})(5,:,3,:))));
    
    %% Logit regression Dummy Segment and X1 (% correct classification)
    
    % define optimization options
    options = optimset('Display','off','FunValCheck','on','MaxFunEvals',Inf, ...
    'MaxIter',Inf,'TolFun',1e-5,'TolX',1e-6);
    
    pct         = zeros(1000,1);
    const       = ones(1250,1); % need a constant regressor too
    warning off
    
    for i = 1:1000
        % define logistic function
        logfun = @(x) 1./(1+exp(-x));
        % construct regressor matrix (data points x regressors)
        y  = SaveSegment.(rc_names{rc_idx})(:,i);
        X1 = SaveX1.(rc_names{rc_idx})(:,i);
        X  = [X1 const];
        % the cost function to minimize is the negative-log-likelihood for the binomial logit
        % here, the eps is a hack to ensure that the log returns a finite number
        costfun     = @(pp) -sum((y .* log(feval(logfun,X*pp')+eps)) + ...
            ((1-y) .* log(1-feval(logfun,X*pp')+eps)));
        x0          = zeros(1,size(X,2));
        params      = fminunc(costfun,x0,options);
        prediction  = X*params'>=0.5;
        % calculate percent correct
        pct(i) = sum(prediction==y) / length(y) * 100;
    end
    warning on
    
    %% Average classification
    AvCorrectClassified     = mean(pct);
    
    %% Coefficients and standard errors
    
    MeanCoefficient         = nanmean(Coefficients.(rc_names{rc_idx}),4);
    EmpiricalStErrors       = nanstd(Coefficients.(rc_names{rc_idx}),0,4);
    MeanStandardErrors      = nanmean(StandardErrors.(rc_names{rc_idx}),4);
    MeanOutsideShare        = mean(OutsideShare.(rc_names{rc_idx}),2);
    CompTime                = sum(ComputationTime.(rc_names{rc_idx}),4);
    
    % Lilliefors Normality test
    % Nesting parameter
    
    [hNL,pNL]                 = lillietest(squeeze(Coefficients.(rc_names{rc_idx})(4,:,2,:))); % Nested logit
    
    % Random coefficient parameter
    
    [hRCL,pRCL]               = lillietest(squeeze(Coefficients.(rc_names{rc_idx})(5,:,3,:))); % RCL
    
    % Jarque-Bera Normality test
    % Nesting parameter
    
    [JBhNL,JBpNL]                 = jbtest(squeeze(Coefficients.(rc_names{rc_idx})(4,:,2,:))); % Nested logit
    
    % Random coefficient parameter
    
    [JBhRCL,JBpRCL]               = jbtest(squeeze(Coefficients.(rc_names{rc_idx})(5,:,3,:))); % RCL
    
    %% AIC and BIC criteria to choose the appropriate model
    
    AICrc   = squeeze(AIC.(rc_names{rc_idx})); % rows=number of estimators; cols=number of MC simulations
    MC      = size(AICrc,2);
    MinAIC  = zeros(1,MC);
    
    for j=1:MC
        [MinAIC(:,j)]     = find(AICrc(:,j)==min(AICrc(:,j)));
    end
    
    BICrc      = squeeze(BIC.(rc_names{rc_idx})); % rows=number of estimators; cols=number of MC simulations
    MinBIC     = zeros(1,MC);
    
    for j=1:MC
        [MinBIC(:,j)]     = find(BICrc(:,j)==min(BICrc(:,j)));
    end
    
    CountAIC1 = countmember(1,MinAIC); % count how many times Logit model is the preferred model according to AIC criterion
    CountBIC1 = countmember(1,MinBIC); % count how many times Logit model is the preferred model according to BIC criterion
    
    CountAIC2 = countmember(2,MinAIC); % count how many times NL model is the preferred model according to AIC criterion
    CountBIC2 = countmember(2,MinBIC); % count how many times NL model is the preferred model according to BIC criterion
    
    CountAIC3 = countmember(3,MinAIC); % count how many times RCL model is the preferred model according to AIC criterion
    CountBIC3 = countmember(3,MinBIC); % count how many times RCL model is the preferred model according to BIC criterion
    
    
    %% Elasticity analysis
    
    MeanOwnElasticity         = nanmean(OwnElast.(rc_names{rc_idx}),4);
    StDevOwnElasticity        = nanstd(OwnElast.(rc_names{rc_idx}),0,4);
    
    MeanCrossElastSameSeg     = nanmean(CrossElastSameSeg.(rc_names{rc_idx}),4);
    StDevCrossElastSameSeg    = nanstd(CrossElastSameSeg.(rc_names{rc_idx}),0,4);
    
    MeanCrossElastDiffSeg     = nanmean(CrossElastDiffSeg.(rc_names{rc_idx}),4);
    StDevCrossElastDiffSeg    = nanstd(CrossElastDiffSeg.(rc_names{rc_idx}),0,4);

    OneOwnElasticity          = nanmean(OneOwnElast.(rc_names{rc_idx}),4);
    OneStDevOwnElasticity     = nanstd(OneOwnElast.(rc_names{rc_idx}),0,4);
    
    OneCrossElasticSameSeg    = nanmean(OneCrossElastSameSeg.(rc_names{rc_idx}),4);
    StDevCrossElasticSameSeg  = nanstd(OneCrossElastSameSeg.(rc_names{rc_idx}),0,4);
    
    OneCrossElasticDiffSeg    = nanmean(OneCrossElastDiffSeg.(rc_names{rc_idx}),4);
    StDevCrossElasticDiffSeg  = nanstd(OneCrossElastDiffSeg.(rc_names{rc_idx}),0,4);
   
    MeanDiversionSameSeg      = nanmean(DiversionSameSeg.(rc_names{rc_idx}),4);
    StDevDiversionSameSeg     = nanstd(DiversionSameSeg.(rc_names{rc_idx}),0,4);
    
    MeanDiversionDiffSeg      = nanmean(DiversionDiffSeg.(rc_names{rc_idx}),4);
    StDevDiversionDiffSeg     = nanstd(DiversionDiffSeg.(rc_names{rc_idx}),0,4);

    MeanOneDiversionSameSeg   = nanmean(OneDiversionSameSeg.(rc_names{rc_idx}),4);
    StDevOneDiversionSameSeg  = nanstd(OneDiversionSameSeg.(rc_names{rc_idx}),0,4);
        
    MeanOneDiversionDiffSeg   = nanmean(OneDiversionDiffSeg.(rc_names{rc_idx}),4);
    StDevOneDiversionDiffSeg  = nanstd(OneDiversionDiffSeg.(rc_names{rc_idx}),0,4);
        
    % Display and save results
    
    variables       = {'const';'dummy';'x1';'pho';'sigma'};
    titleCoeff      = {'Coefficients'};
    titleEstimator  = {'True','Logit','NL','RCL'};
    resultsCoeff    = [Coefftrue.(rc_names{rc_idx}), squeeze(MeanCoefficient)];
    
    disp(titleCoeff);
    disp([variables,num2cell(resultsCoeff)])
    
    variables       = {'const';'dummy';'x1';'pho';'sigma'};
    titleStErr1     = {'Empirical StError'};
    resultsStErr1   = squeeze(EmpiricalStErrors);
    
    disp(titleStErr1);
    disp([variables,num2cell(resultsStErr1)])
    
    variables       = {'const';'dummy';'x1';'pho';'sigma'};
    titleStErr2     = {'Calculated StError'};
    resultsStErr2   = squeeze(MeanStandardErrors);
    
    disp(titleStErr2);
    disp([variables,num2cell(resultsStErr2)])
    
    xlswrite2007(File,titleCoeff, [rc_names{rc_idx},TitleSheet],'A1')
    xlswrite2007(File,titleEstimator,[rc_names{rc_idx},TitleSheet],'B1')
    xlswrite2007(File,resultsCoeff,[rc_names{rc_idx},TitleSheet],'B2')
    xlswrite2007(File,variables,[rc_names{rc_idx},TitleSheet],'A2')
    
    xlswrite2007(File,titleStErr1,[rc_names{rc_idx},TitleSheet],'A8')
    xlswrite2007(File,titleEstimator,[rc_names{rc_idx},TitleSheet],'B8')
    xlswrite2007(File,resultsStErr1,[rc_names{rc_idx},TitleSheet],'C9')
    xlswrite2007(File,variables,[rc_names{rc_idx},TitleSheet],'A9')
    
    xlswrite2007(File,titleStErr2,[rc_names{rc_idx},TitleSheet],'A15')
    xlswrite2007(File,titleEstimator,[rc_names{rc_idx},TitleSheet],'B15')
    xlswrite2007(File,resultsStErr2,[rc_names{rc_idx},TitleSheet],'C16')
    xlswrite2007(File,variables,[rc_names{rc_idx},TitleSheet],'A16')
    
    variables  = {'mins0';'avgs0';'maxs0'};
    xlswrite2007(File,variables,[rc_names{rc_idx},TitleSheet],'A23')
    xlswrite2007(File,squeeze(MeanOutsideShare),[rc_names{rc_idx},TitleSheet],'B23')
    
    variables  = {'CPUTime'};
    xlswrite2007(File,variables,[rc_names{rc_idx},TitleSheet],'A27')
    xlswrite2007(File,(squeeze(CompTime))',[rc_names{rc_idx},TitleSheet],'C27')
    
    % Elasticities
    %% Mean values
    
    variables  = {'Segment1';'Segment2'};
    titleElast = {'Own Elasticity'};
    
    xlswrite2007(File,titleElast,[rc_names{rc_idx},TitleSheet],'A29')
    xlswrite2007(File,variables,[rc_names{rc_idx},TitleSheet],'A30')
    xlswrite2007(File,squeeze(MeanOwnElasticity(:,1,:)),[rc_names{rc_idx},TitleSheet],'C30')
    
    titleElast = {'Cross Elasticities Same Seg'};
    
    xlswrite2007(File,titleElast,[rc_names{rc_idx},TitleSheet],'A33')
    xlswrite2007(File,variables,[rc_names{rc_idx},TitleSheet],'A34')
    xlswrite2007(File,squeeze(MeanCrossElastSameSeg(:,1,:)),[rc_names{rc_idx},TitleSheet],'C34')
    
    titleElast = {'Cross Elasticities Diff Seg'};
    
    xlswrite2007(File,titleElast,[rc_names{rc_idx},TitleSheet],'A37')
    xlswrite2007(File,variables,[rc_names{rc_idx},TitleSheet],'A38')
    xlswrite2007(File,squeeze(MeanCrossElastDiffSeg(:,1,:)),[rc_names{rc_idx},TitleSheet],'C38')
    
    %% AIC and BIC
    
    titleAIC    = {'AIC'};
    xlswrite2007(File,titleAIC,[num2str((rc_names{rc_idx})),TitleSheet],'A41')
    xlswrite2007(File,titleEstimator,[num2str((rc_names{rc_idx})),TitleSheet],'B41')
    xlswrite2007(File,[CountAIC1,CountAIC2,CountAIC3],[num2str((rc_names{rc_idx})),TitleSheet],'C42')
    
    titleBIC    = {'BIC'};
    xlswrite2007(File,titleBIC,[num2str((rc_names{rc_idx})),TitleSheet],'A44')
    xlswrite2007(File,titleEstimator,[num2str((rc_names{rc_idx})),TitleSheet],'B44')
    xlswrite2007(File,[CountBIC1,CountBIC2,CountBIC3],[num2str((rc_names{rc_idx})),TitleSheet],'C45')
    
    %% Lilliefors Normality tests
    
    title    = {'Lilliefors Normality test Nest'};
    xlswrite2007(File,title,[num2str((rc_names{rc_idx})),TitleSheet],'A47')
    xlswrite2007(File,titleEstimator,[num2str((rc_names{rc_idx})),TitleSheet],'B47')
    xlswrite2007(File,[NaN,hNL,NaN;NaN,pNL,NaN],[num2str((rc_names{rc_idx})),TitleSheet],'C48')
    
    title    = {'Lilliefors Normality test RC'};
    xlswrite2007(File,title,[num2str((rc_names{rc_idx})),TitleSheet],'A51')
    xlswrite2007(File,titleEstimator,[num2str((rc_names{rc_idx})),TitleSheet],'B51')
    xlswrite2007(File,[NaN,NaN,hRCL;NaN,NaN,pRCL],[num2str((rc_names{rc_idx})),TitleSheet],'C52')
    
    %% Jarque-Bera Normality tests
    
    title    = {'Jarque-Bera Normality test Nest'};
    xlswrite2007(File,title,[num2str((rc_names{rc_idx})),TitleSheet],'A54')
    xlswrite2007(File,titleEstimator,[num2str((rc_names{rc_idx})),TitleSheet],'B54')
    xlswrite2007(File,[NaN,JBhNL,NaN;NaN,JBpNL,NaN],[num2str((rc_names{rc_idx})),TitleSheet],'C55')
    
    title    = {'Jarque-Bera Normality test RC'};
    xlswrite2007(File,title,[num2str((rc_names{rc_idx})),TitleSheet],'A58')
    xlswrite2007(File,titleEstimator,[num2str((rc_names{rc_idx})),TitleSheet],'B58')
    xlswrite2007(File,[NaN,NaN,JBhRCL;NaN,NaN,JBpRCL],[num2str((rc_names{rc_idx})),TitleSheet],'C59')
    
    %% Diversion Ratios
    
    variables  = {'Segment1';'Segment2'};
    titleElast = {'Diversion Ratio Same Segment'};
    
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A63')
    xlswrite2007(File,variables,[num2str((rc_names{rc_idx})),TitleSheet],'A64')
    xlswrite2007(File,squeeze(MeanDiversionSameSeg(:,1,:)),[num2str((rc_names{rc_idx})),TitleSheet],'C64')
    
    titleElast = {'Diversion Ratio Diff Segment'};
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A67')
    xlswrite2007(File,variables,[num2str((rc_names{rc_idx})),TitleSheet],'A68')
    xlswrite2007(File,squeeze(MeanDiversionDiffSeg(:,1,:)),[num2str((rc_names{rc_idx})),TitleSheet],'C68')
    
    titleElast = {'One Diversion Ratio Same Segment'};
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A71')
    xlswrite2007(File,(squeeze(MeanOneDiversionSameSeg(:,1,:)))',[num2str((rc_names{rc_idx})),TitleSheet],'C72')
    titleElast = {'One Diversion Ratio Diff Segment'};
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A74')
    xlswrite2007(File,(squeeze(MeanOneDiversionDiffSeg(:,1,:)))',[num2str((rc_names{rc_idx})),TitleSheet],'C75')
        
    %% correctly classied
    xlswrite2007(File,{'correct classified'},[num2str((rc_names{rc_idx})),TitleSheet],'A77')
    xlswrite2007(File,AvCorrectClassified,[num2str((rc_names{rc_idx})),TitleSheet],'C77')
    
    %% Standard deviation elasticities
    
    variables  = {'Segment1';'Segment2'};
    titleElast = {'Standard deviation Own Elasticity'};
    
    xlswrite2007(File,titleElast,[rc_names{rc_idx},TitleSheet],'A79')
    xlswrite2007(File,variables,[rc_names{rc_idx},TitleSheet],'A80')
    xlswrite2007(File,squeeze(StDevOwnElasticity(:,1,:)),[rc_names{rc_idx},TitleSheet],'C80')
    
    titleElast = {'Standard deviation Cross Elasticities Same Seg'};
    
    xlswrite2007(File,titleElast,[rc_names{rc_idx},TitleSheet],'A83')
    xlswrite2007(File,variables,[rc_names{rc_idx},TitleSheet],'A84')
    xlswrite2007(File,squeeze(StDevCrossElastSameSeg(:,1,:)),[rc_names{rc_idx},TitleSheet],'C84')
    
    titleElast = {'Standard deviation Cross Elasticities Diff Seg'};
    
    xlswrite2007(File,titleElast,[rc_names{rc_idx},TitleSheet],'A87')
    xlswrite2007(File,variables,[rc_names{rc_idx},TitleSheet],'A88')
    xlswrite2007(File,squeeze(StDevCrossElastDiffSeg(:,1,:)),[rc_names{rc_idx},TitleSheet],'C88')
    
   %% Standard deviation diversion ratio
    
    variables  = {'Segment1';'Segment2'};
    titleElast = {'St.Dev. Diversion Same Segment'};
    
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A91')
    xlswrite2007(File,variables,[num2str((rc_names{rc_idx})),TitleSheet],'A92')
    xlswrite2007(File,squeeze(StDevDiversionSameSeg(:,1,:)),[rc_names{rc_idx},TitleSheet],'C92')
    
    titleElast = {'Standard deviation Diversion Ratio Diff Segment'};
    
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A95')
    xlswrite2007(File,variables,[num2str((rc_names{rc_idx})),TitleSheet],'A96')
    xlswrite2007(File,squeeze(StDevDiversionDiffSeg(:,1,:)),[rc_names{rc_idx},TitleSheet],'C96')
    
    titleElast = {'Standard deviation One Diversion Ratio Same Segment'};
    
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A99')
    xlswrite2007(File,(squeeze(StDevOneDiversionSameSeg(:,1,:)))',[rc_names{rc_idx},TitleSheet],'C100')
    
    titleElast = {'Standard deviation One Diversion Ratio Diff Segment'};
    
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A103')
    xlswrite2007(File,(squeeze(StDevOneDiversionDiffSeg(:,1,:)))',[rc_names{rc_idx},TitleSheet],'C104')   
    
    % One elasticities
   
    titleElast = {'One Own Elasticity'};
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A105')
    xlswrite2007(File,(squeeze(OneOwnElasticity(:,1,:)))',[rc_names{rc_idx},TitleSheet],'C106')
    
    titleElast = {'One Cross Elasticity Same Segment'};
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A108')
    xlswrite2007(File,(squeeze(OneCrossElasticSameSeg(:,1,:)))',[rc_names{rc_idx},TitleSheet],'C109')

    titleElast = {'One Cross Elasticity Diff Segment'};
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A111')
    xlswrite2007(File,(squeeze(OneCrossElasticDiffSeg(:,1,:)))',[rc_names{rc_idx},TitleSheet],'C112')    
   
    % One standard deviations elasticities
    
    titleElast = {'Standard deviation One Own Elasticity'};
    
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A116')
    xlswrite2007(File,(squeeze(OneStDevOwnElasticity(:,1,:)))',[rc_names{rc_idx},TitleSheet],'C117')
    
    titleElast = {'Standard deviation One Cross Elasticity Same Seg'};
    
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A120')
    xlswrite2007(File,(squeeze(StDevCrossElasticSameSeg(:,1,:)))',[rc_names{rc_idx},TitleSheet],'C121')
    
    titleElast = {'Standard deviation One Cross Elasticity Diff Seg'};
    
    xlswrite2007(File,titleElast,[num2str((rc_names{rc_idx})),TitleSheet],'A123')
    xlswrite2007(File,(squeeze(StDevCrossElasticDiffSeg(:,1,:)))',[rc_names{rc_idx},TitleSheet],'C124')      
end

ExcelWorkbook.Save
ExcelWorkbook.Close(false) % Close Excel workbook.
Excel.Quit;
delete(Excel);